Spherical Collapse and Cluster Counts in Modified Gravity Models 



Matthew C. Martino, Hans F. Stabenau & Ravi K. ShethEl 
Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19103 

Modifications to the gravitational potential affect the nonlinear gravitational evolution of large 
scale structures in the Universe. To illustrate some generic features of such changes, we study the 
evolution of spherically symmetric perturbations when the modification is of Yukawa type; this is 
non-trivial, because we should not and do not assume that Birkhoff 's theorem applies. We then show 
how to estimate the abundance of virialized objects in such models. Comparison with numerical 
simulations shows reasonable agreement: When normalized to have the same fluctuations at early 
times, weaker large scale gravity produces fewer massive halos. However, the opposite can be true for 
models that are normalized to have the same linear theory power spectrum today, so the abundance 
of rich clusters potentially places interesting constraints on such models. Our analysis also indicates 
that the formation histories and abundances of sufficiently low mass objects are unchanged from 
standard gravity. This explains why simulations have found that the nonlinear power-spectrum at 
large k is unaffected by such modifications to the gravitational potential. In addition, the most 
massive objects in CMB-normalized models with weaker gravity are expected to be similar to the 
high-redshift progenitors of the most massive objects in models with stronger gravity. Thus, the 
difference between the cluster and field galaxy populations is expected to be larger in models with 
stronger large-scale gravity. 



I. INTRODUCTION 

When we consider the cosmological data from WMAP, 
supernovae la, galaxy clustering on large scales, and 
cross-correlations between galaxies and the CMB, we are 
faced with several possible conclusions. One is that we 
live in a spatially flat Friedmann-Robertson- Walker uni- 
verse currently dominated by either a cosmological con- 
stant or repulsive dark energy. The best fit for the di- 
mensionless energy density parameters are VL m = 0.28 
and J7a = 0.72, within the concordance ACDM Model 
[l| . The advantage of this model is that nearly all avail- 
able observational data supports it; the disadvantage is 
that it requires the vast majority of the energy density 
in the universe to be in two unknown substances, dark 
matter and dark energy. 

This conclusion rests on the accuracy of our current 
gravity model, General Relativity. The key equation of 
General Relativity, Einstein's equation, relates the curva- 
ture and the expansion rate of the Universe to its matter 
and energy content. The current paradigm is to mod- 
ify the matter content of the universe, by including dark 
matter and dark energy, to account for observations. In- 
stead, however, we might modify how the universe curves 
in response to matter, which would mean modifying our 
theory of gravity. 

There are many possible ways to modify gravity, de- 
pending on what one wishes to "fix" . For example, 
MOND @,0| removes the need for dark matter to account 
for galactic rotation curves and has several other interest- 
ing results, but seems to fail at the scale of galaxy clus- 
ters, with even optimistic accountings needing roughly as 
much dark matter as baryonic matter. At the other end 



is something like DGP or conformal gravity, which hopes 
to account for the acceleration of the universe without 
invoking a cosmological constant 0, H, H, 0] • In addition, 
there are other models which seek to unify dark matter 
and dark energy [1, [1, EE EI El • What we seek to do 
in this paper is to study the problem more phenomolog- 
ically. For example, regardless of how the force law for 
gravity is modified, it will often be stronger or weaker, 
relative to the standard model, at larger or shorter scales. 
One way to parameterize this is to introduce a modified 
Yukawa-like potential for a point mass: 



(r) = Gm 



I + a(I-e- r / r =) 



(1) 
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EH EE EB| where a indicates the strength and r s 
the scale (constant in physical rather than comoving co- 
ordinates) on which this modification is most relevant. 
On scales smaller than r s , 0(r) reduces to the standard 
Newtonian potential; on larger scales it transitions to the 
Newtonian potential multiplied by a factor of (1 + a). 
This is similar to the interaction considered in [l7j |. in 
which a long range dark matter interaction is introduced 
yielding a different Yukawa-like potential that instead 
modifies gravity on short length scales. 

Note that this is not a cosmological model: there is 
no prescription for determining things like the expansion 
factor and the resulting Hubble factor. Like previous 
workers, we will assume that these are the same as in 
the standard cosmology. The main goal of studying such 
a model is to gain intuition for some generic effects of 
modifications to standard gravity. 

For example, in the linear theory of such a model, the 
growth of fluctuations is fc-dependent - it is not in 
standard gravity. As a result, a smooth spherical region 
within which the density is the same as the background 
universe will evolve. This qualitatively different behavior 
fcdam standard gravity has not been emphasized - so it 
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is worth showing the argument explicitly. Consider the 
density field smoothed on scale R at some early time ti. 
We can write this field in terms of its Fourier modes, and 
the (Fourier Transform of the) smoothing kernel as 



S R (x, ti) = J dk exp(ik • x) S(k) W(kR). 



(2) 



The linearly evolved field is 



<5fl(x,i) 



J e?k exp(ik • x) 



D(k, t) 
D(k,U) 



S(k)W(kR), (3) 



where D is the linear theory growth factor. In stan- 
dard gravity, D is independent of k, so if <5r(x, fcj) = 
then <5r(x, t) = also. But if D depends on k, then if 
<5r(x) = at some time t, it will, in general, be non-zero 
at other times (the exception being if the fc-dependence 
of W happens to exactly cancels that of D). For the 
Yukawa-like modification considered here, the k depen- 
dence of the spherical top hat filter does not cancel that of 
D. Thus, we are led to the rather remarkable conclusion 
that, when the gravitational potential has been modi- 
fied, then linear theory predicts that a spherical tophat 
patch within which the density is the same as the back- 
ground will evolve! The reason why can be traced to the 
fact that Birkhoff 's Theorem no longer applies once the 
Newtonian potential has been modified. Without this 
Theorem, the spherical top hat filter is no longer special, 
and our common sense prejudice from standard gravity 
- that initially overdense regions become denser, under- 
dense regions less dense, but regions within which the 
density is the same as the background do not evolve - 
must be treated with caution. 

The evolution of nonlinear clustering in this model has 
been studied using numerical simulations by [THEH . Our 
goal in what follows is to provide a framework for under- 
standing this nonlinear evolution more completely. To 
this end, we will use the spherical evolution model, which 
has found extensive use in the standard model - it is used 
to motivate estimates of the abundance of nonlinear ob- 
jects [l8j , a crucial ingredient in methods which describe 
the growth of nonlinear gravitational clustering [19(. It 
also provides a framework for discussing the environmen- 
tal dependence of clustering [2(| [2l| . 

Section |TT] summarizes what is known from the linear 
theory of this Yukawa-like model, and then shows our 
estimate of the key, and sometimes subtle, changes to the 
spherical evolution model. Scction[TTT]compares this work 
to the simulations of [la. A final section summarizes. 



II. THE MODEL 

A. Linear theory and a slightly modified potential 

We begin by considering the evolution of density per- 
turbations. This can be done by either considering the 



fluid equations in expanding coordinates, or by consid- 
ering the conservation of stress-energy V a T ab = 0. If 
we start with a smooth background, add small pertur- 
bations, and linearize the resulting equations, we get a 
second order differential equation for the evolution of the 
fractional overdensity, 8(x, t) that depends on time, scale 
factor a, the Hubble parameter H = a/a, and the po- 
tential <fi. In standard gravity, we would use the Poisson 
equation to set the relationship between <fi and 5, but here 
we will assume a modified Poisson equation that results 
in the above potential, equation |T|). 



6 + 2H(t)S = V 2 



(4) 



It is easier to work with the Fourier transform of this 
equation: 



2HS k = - 



1 



1 - 



H 2 



H 2 S k 



(5) 

This can be solved relatively easily to determine the lin- 
ear growth of structure associated with equation (fTJ) . 

The dashed lines in Figure [1] show how this growth dif- 
fers from that in the standard model. Note in particular 
that S^ in (t) = D(k, t) <5[ nitial , whereas the standard model 
has no k dependence in the growth factor D(t). The fig- 
ure shows the effects that one expects to see, namely 
that for negative a the growth factor is smaller at small 
k (large scales), whereas for positive a the growth factor 
is larger at large scales. Both have a region where they 
deviate from being scale independent, until on very large 
scales they return to scale independence, though with a 
different value than in General Relativity. 

While this is not a significant problem for linear the- 
ory, we decided to explicitly force the potential back to 
normal gravity at large scales for reasons which will be- 
come clear in the next section, briefly, because we want 
to assume that the cosmological model is indistinguish- 
able from ACDM, on the largest scales we want gravity 
to be the same as in ACDM. The potential we are using 
from this point onward is: 



Gm 



l + a(l - e- r / r =) - a(l 



-r/r c 



(6) 



With this potential, the r c 3> r s term serves to make 
explicit the return to normal gravity at large scales. The 
linear theory equation becomes 



2HS k 



1 



fc 2 r 



1 - 



H 2 



(7) 



and the linear growth associated with this solution 
is shown as the solid lines in Figure [TJ when r c — 
IQh^ 1 Mpc. The dotted line shows what happens if 
r c = 350/i _1 Mpc - a scale which is large compared to 
that probed by BAOs. Although the analysis which fol- 
lows uses r c = 70/i _1 Mpc, the results which follow are 
not sensitive to this choice. 



3 





1 1 1 1 1 


10 










y\ ] 

/ / \ 

/ / \ 


1 





Id" 6 


10" s 0.0001 0.001 0.01 0.1 
k (h/Mpc) 


1 




Sb.i 

a 

\ 




\ / 

\ / ] 


0.01 




1 






io-° 


10" 5 0.0001 0.001 0.01 0.1 



k (h/Mpc) 

FIG. 1: Ratio of linear theory growth factor to that in stan- 
dard ACDM, at a = 1, when r B = 5/i _1 Mpc and a = 1 (top) 
and a = — 1 (bottom). Dashed lines show this ratio for the 
model in equation (fTJ) , and solid lines for equation ((6)1 . For the 
solid lines, r c = 7Qh~ x Mpc, and dotted r c — 350/i." 1 Mpc. 



B. Spherical collapse 

Excursion set methods [13, HE H3 are used to esti- 
mate the number density of collapsed haloes, the merger 
rates of haloes, the conditional mass function of progen- 
itors as a function of final halo mass and time, and the 
nonlinear counts-in-cells distribution, all of which help us 
to link what we observe about the properties of galaxies 
and galaxy clusters with cosmology. Essentially, excur- 
sion set methods relate the properties of halos today to 
the initial density fluctuation field. An advantage of such 
methods is that we only need a few things in order to use 
them: one is an assumption that the initial fluctuations 
are small and Gaussian, the other is a model for deter- 
mining how dense something must have been initially 
to collapse at a given time. The first is given to us by 
WMAP, the second is more difficult. 

A simple model for how to determine this critical den- 
sity is given by the spherical evolution model, in which 
one considers a spherical tophat perturbation in the ini- 



tial density fluctuation field. In the standard model, the 
gravitational evolution of such a patch is determined only 
by the mass within it, and so one can determine how 
overdense such a patch needs to be initially in order to 
collapse by a given time. This critical overdensity S c 
generically depends on the background cosmology [251 ]. 

The spherical collapse calculation begins with the 
statement that the force driving the acceleration is re- 
lated to the gravitational potential by 

^ = F=-V$(r). (8) 
This can be integrated once to get 

where <&(r) is the integral of the potential over the mass 
distribution, and C is the total energy of the patch, which 
is constant. In standard gravity, the potential of a shell 
of mass M is the same as that of a point mass at the 
center of the sphere, so $(r) reduces to GM(< r)/r. 
The constant C can be related to the initial overdensity 
and/or expansion rate of the patch: the initial expan- 
sion rate is given by the Hubble expansion rate of the 
background in which the patch is embedded, namely in 
comoving coordinates ii = 0, so fi — diXi = di/a,i(aiXi), 
so (dr/dt)i = HiVi. Including a cosmological constant 
presents no conceptual difference. 

In standard gravity one can directly solve this equa- 
tion. The solution is a cycloid for which the critical value 
of the initial overdensity required for collapse, S c , docs 
not depend on the initial size of the patch. This scale 
independence of 5 C is a result of Birkhoff 's Theorem: the 
evolution of a tophat sphere is the same as that given 
by Friedmann's equations, so the actual size of the patch 
drops out. 

When gravity is modified, things are no longer so sim- 
ple. For example, when the potential is given by equa- 
tion ©, Birkhoff 's theorem no longer applies: A parti- 
cle offcenter in a uniform spherical shell will feel a force 
from the shell because we no longer have a 1 /r 2 force law. 
This has two consequences. First, equation (J8]) can still 
be integrated once to get equation ©, only now $(r) has 
contributions from both the internal and external mass 
distributions. We can get $ by integrating equation ^ of 
the mass distribution, leaving C and the initial value for 
dr/dt to be determined. As before, C is the total energy 
(constant in time), and we set (dr/dt)i = HiTi. (This 
was why we used equation [5] rather than equation [TJ) 
Second, whereas evenly spaced concentric shells remain 
evenly spaced in the standard tophat model, this is no 
longer the case when the potential is modified. As a re- 
sult, the initial tophat perturbation develops a nontrivial 
density profile as it collapses. 

Of course, neither of these changes prevent us from 
solving for the evolution of r(t). Our main interest in 
what follows is not in the details of how the density profile 
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FIG. 2: Convergence to solution as number of shells increases 
for two objects with mass 1O 14 - 5 /i _1 M , with a = 0.5 (top) 
and a — 1.0 (bottom). The thick solid lines are actually the 
positions of the shell on the edge of the density perturbation 
for 3,10, 20, 40, 50, 60, 80, 100, and 200 shells, note they all 
lie nearly on top of each other. The dashed curves show the 
evolution when a — 0. 
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FIG. 3: Ratio of initial density required for collapse at the 
present time to that in the standard gravity model, when 
the background cosmology is ACDM, r s = 5/i _1 Mpc and 
r c = 70/i _1 Mpc. From top to bottom, curves show models in 
which q = — 1, —0.5, 0, 0.5 and 1 (note that a = is standard 
gravity). 



to account for the lack of Birkhoff 's theorem, regardless 
of a or r s . We also found that for objects of mass up 
to 1O 1 5/i _1 M0, using 3 (linearly spaced) wide shells pro- 
vided a 5 C that was within 1% of 200 shells (see Figure^. 
For higher mass objects we found that we needed more 
than 3 shells, but by 40 shells S c is within 0.02% of S c 
calculated with 200 shells. 



is modified (this is interesting in its own right), but in the 
modification to the critical density required for collapse. 
To estimate this in practice, two things require care, both 
of which are related to the breakdown of Birkhoff 's the- 
orem. The first is that, because the shape of the per- 
turbation evolves, one must follow the evolution not just 
of a shell at the boundary of the perturbation, but of 
a series of concentric shells. So the question is: How 
finely spaced must the shells be before one converges to 
the correct solution? The second is that one now cares 
not only about the mass initially interior to the initial 
boundary, but the mass exterior as well. In this case, the 
question is: How far beyond the initial boundary of the 
perturbation matters before one reaches convergence? 

Therefore we start with an initial patch which is sub- 
stantially larger than that within which there is an initial 
overdensity, and use a simple 1-dimensional N-body sim- 
ulation to solve for r(t). We found that volumes having 
twice the initial comoving radius or larger were sufficient 



Having determined that our numerical solution had 
converged, we evaluated r(t) at the present time in an 
ACDM background model, for a grid of initial sizes and 
overdensities. By finding which pairs of initial density 
and size Ti when evolved result in r (to) = 0, we obtained 
S c i{fi)- Since the initial overdensity is always small, we 
can use the fact that M = (4n/S)rfpi(l + S t ) w (47r/3)rf 
to express this critical density as a function of mass rather 
than initial radius. This is shown in Figure [3l Notice 
that 5 C depends on mass; this is not unexpected, be- 
cause patches which remain smaller than r s throughout 
their evolution (and become small mass halos) are un- 
likely to notice any modification, whereas those which 
are larger than r s at any time during their evolution will. 
This mass-dependence of S c means that we expect to see 
a variation in cluster abundances only at masses larger 
than ~ 10 14 /i _1 Mq. This is a consequence of the fact 
that r s = hh' 1 Mpc, for which the mass scale is about 
5 x lO 13 /^ 1 M . (For fixed a, the mass scale is propor- 
tional to r 3 .) 
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FIG. 4: (Color online) Expected halo abundance, 
log dN/d log m, as a function of mass m for models with a = 1 
(short dashed) and a = — 1 (long dashed), the solid line is for 
a = 0. The rapid rise of the negative a barrier in figure [3] 
is why there is a larger effect on the abundance of high mass 
halos for negative a than positive a. 

C. The abundance of virialized objects 

Figure [3] shows that, when the potential is modified, 
then S c is no longer scale-independent. Because it de- 
pends on mass, the relevant excursion set problem is one 
with a 'moving' rather than 'constant' barrier, so it is 
of the type first studied by [2l|. Figure H] shows the re- 
sult of using this formalism to estimate the abundance of 
virialized objects. Briefly, making this estimate requires 
that one generate an ensemble of random walks in the 
(Si, Si) plane, where Si = of(Af) is the variance in the 
initial fluctuation field when smoothed on scale T{. Since 
CTi(M) is a monotonic function of M, the variables Si, 
M and r, are essentially equivalent to one another. In 
particular, 

*-/f i ^ W (10) 

where W(x) = (3/x 3 ) (since — xcosx). One then 
finds the first crossing distribution f(Si)dSi of the 'bar- 
rier' 5 C i(M) = S C i(Si). The abundance of objects is 
dn/dlnMdlnM = (p/M) f(S i )dS l . 

Figure H] shows that the abundance of massive ob- 
jects increases as a increases, whereas the abundance 
of low mass objects is essentially unchanged from when 
a = 0. We argued above that this makes intuitive sense: 
the lower mass objects do not feel the change in grav- 
ity because they were smaller than r s throughout their 
evolution; the more massive objects are able to become 



even more massive if a is positive, since then gravity is 
stronger. 



D. Choice of normalization and incompatibility 
with standard gravity 

Before moving on, it is worth noting that we were care- 
ful to describe and implement the excursion set approach 
in initial rather than linearly evolved variables. In stan- 
dard gravity (a = 0), it is more common to phrase the 
discussion in terms of linearly evolved variables. Since 
the linear growth factor is independent of k when a = 0, 
this is straightforward. However, this is no longer the 
case when a/0, because of the fc-dependence in D(k, t). 

Indeed, the barrier shape in Figure [3] is qualitatively 
like that of the linear growth factor in Figure [TJ so one 
might ask if the difference that we see in the mass func- 
tion is entirely a consequence of the fc-dependence of the 
linear growth factor. More specifically, the differences 
shown in Figure [4] are really a consequence of two effects: 
first, we have assumed that Si(M) is the same for all 
a. Therefore, So(M) is not: the linear theory evolution 
when a > results in more large scale power than when 
ot < 0, so the rms fluctuation in the present day fields 
are different - in the jargon, erg at z — is larger for the 
a > models. Since we know that, in standard gravity, 
the abundance of massive halos is exponentially sensitive 
to us j one might wonder if this alone accounts for much of 
the effect. (Later on, we will discuss another consequence 
of normalizing the models at the initial rather than final 
time.) To quantify this, we would like to compare the 
predicted abundances when the models are normalized 
so that Sq(M), rather than Si(M), is the same. This 
will isolate the effect of a mass-dependent S c (M\a, r s , r c ) 
on the halo abundances. 

Therefore, we used the excursion set approach with 
the same constant barrier height as one would have in 
the standard (a — 0) linearly evolved gravity model, and 
then evaluated Sq(M), rather than Si(M), using the k- 
dependent linear growth factor. I.e., 

So^| fD^k,t )^^wHkn), (11) 

The resulting first crossing distribution f(So)dSo is the 
same as in standard gravity (after all the barrier is con- 
stant), but when expressed as a function of mass M, the 
abundances are different because the relation between So 
and M depends on (a, r s , r c ) . The dashed line in Figured 
shows that this method yields a mass function that also 
differs substantially from the standard one: it drops sub- 
stantially below unity for large M . In this case, however, 
the drop is entirely due to the fc-dependent growth fac- 
tor. In addition, notice that although it is qualitatively 
like the solid line, for which 5 C is mass dependent, it can 
be substantially different (i.e., the ratio of the solid to 
the dashed line is greater than unity) at large M. This 
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FIG. 5: Ratio of halo abundances when a = — 1 (top) and 
a = 0.5 (bottom) to that when a — 0. Solid curve shows this 
ratio when the excursion set method is correctly implemented, 
using the initial fluctuation field values Si(M) and a moving 
barrier S C (M), and Si(M) is independent of a. Dashed curve 
uses 5*o (M) from the linearly evolved field and a constant 
barrier S c — 1.686. Dotted curve uses Si(M) and <5 C (M), but 
now Si(M) is modified so that So(M) is the same for both 
values of a. Note that all these approaches produce different 
results. 
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FIG. 6: (Color online) Halo abundances with r s = 5 Mpc 
ratioed to a — 0. In all panels, data points are from simula- 
tion, whereas curves are from our excursion set calculation. In 
the top left panel, (blue) long dashed and (red) short dashed 
curves are for a = 0.5 and —1, respectively, and all models 
had the same initial fluctuation spectrum (so that a = has 
as = 1 today). The top right panel shows the ratio of counts 
in two a = runs, but with different initial conditions, tay- 
lored so that as at z = is the same as in the a = 0.5 (long 
dashed) and —1 (short dashed) runs. Short dashed curves in 
bottom panel show the ratio of the a = 0.5 counts to that in 
the a = run when it has been normalized to have the same 
(linear theory) power spectrum at z — as the a = 0.5 run. 
Long dashed curve shows a similar analysis of the a = — 1 
case. This panel shows the ratio of the numerators in the 
previous two panels. 



shows that there is more to the change in the mass func- 
tion than simply the change in the relation between S 
and mass. 

Whereas the solid line shows results in which the ini- 
tial power spectra are the same for all a (i.e., Si(M) 
is the same for all models), the dashed line shows what 
happens if we adjust the shape of the initial power spec- 
trum in the a — model so that Sq(M) is the same as 
for the a = — 1 model. The dotted line shows the re- 
sult of adjusting instead the initial P(k) of the a = — 1 
model so that it produces the same So(M) as the orig- 
inal a = calculation. In this case, the a = — 1 ini- 
tial conditions now have substantially more large scale 
power, so the predicted abundance of massive halos is 
larger, until the mass dependence of 5 C begins to matter 
(this is not evident in the figure, because it happens at 



M > 10 16 - 7 M Q /h). 

The fact that none of the curves shown in Figure [S] are 
unity for all M, nor are any two curves the same, means 
that cluster abundances in modified gravity models can- 
not be mimicked in standard gravity simply by changing 
the shape of the initial power spectrum so that it agrees 
with modified linear theory at z = 0. Trading 'CMB'- 
normalization for 'cluster' normalization does not work, 
because the cluster mass function depends on the nonlin- 
ear physics of gravitational collapse: Cluster counts are 
sensitive to more than the change to linear theory. For 
CMB-normalized models, we feel that the appropriate 
estimate of the effect of modifying gravity is shown by 
the solid line. In the following section we use numerical 
simulations to test this prediction. 
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III. COMPARISON TO SIMULATIONS 

We now compare our spherical collapse predictions for 
halo abundances with measurements in the simulations 
of [l6(. These simulations followed the evolution of 128 3 
particles in a periodic box of size 100/i _1 Mpc, for various 
choices of a and r s . In all cases, the background cosmol- 
ogy was flat ACDM with fl = 0.3, and the particle mass 
was 1.1 x 10 10 Mq. In addition, the simulations were al- 
ways started from the same initial phases, a feature we 
will exploit shortly. We identify halos in the standard 
way using a friends-of-friends algorithm with link length 
0.2 times the interparticle separation. In what follows, 
we show results from the r s = 5/i _1 Mpc runs. The 
a = simulation, with standard initial conditions has 
ag = 1.0 at z = 0, the corresponding runs for a = 0.5 
and a = — 1 have as = 1.2 and 0.7 respectively. Follow- 
ing our discussion of how the counts depend on the shape 
and normalization of the initial power spectrum, we also 
study results from a — simulations in which the initial 
power spectrum was modified so that, at z = 0, it has 
the same shape as the two a / cases (as = 1.2 and 
0.7). 

Figure [5] shows how the mass function depends on a. 
The panels shows the ratio to a = 0, and curves show 
the predictions from our excursion set with moving or 
standard barrier calculation with standard or modified 
initial power spectrum. The calculation is in reasonably 
good agreement - note in particular that it captures the 
sense of the trends with a. 

Because the simulations all had the same initial phases, 
we were able to perform a slightly more stringent test. 
Namely, we directly compared the masses of individual 
halos in the a — 0.5 and a = — 1 runs with those when 
a = (i.e. standard gravity). The hllcd triangles in 
Figure [7] show the result of selecting the most massive 
halos in the a ^ runs and plotting the number of par- 
ticles they contain versus the number of particles they 
contained in the a = run. The open squares show the 
result of making the selection in the a = run. The 
top and bottom panels show results for a = — 1 and 0.5 
respectively. This illustrates that with stronger gravity 
(larger a), a given halo is more likely to become more 
massive, but this only matters for halos more massive 
than about pr^ . 

IV. CONCLUSIONS 

We presented a study of nonlinear effects in a model 
with a modified gravitational potential (equation [6]) . In 
particular, we showed how the spherical evolution model 
is modified, and the effect this has on the abundance 
of virialized objects. Halos are more massive in mod- 
els where gravity is stronger on large scales (Figure [7]) , 
although this effect is only important for sufficiently mas- 
sive objects whose evolution brings them close to the 
scale r s on which gravity was modified. The effect this 
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FIG. 7: The panel on the top is for a — —1.0 vs. standard, 
the bottom is a = 0.5 vs. standard. The points come from 
selecting big haloes in each realization of the modified gravity 
simulation and then finding the corresponding halo in the 
standard gravity simulation, and then the reverse. N is the 
number of particles in the identified halos. The tilt shows 
that in stronger gravity massive halos are likely to have more 
particles. The "tails" in the lower left are unimportant. 



has on the abundance of massive objects depends on 
how the models are normalized. If normalized so that 
the fluctuation field is the same at early-times (CMB- 
normalized), then the models with a > have more 
massive halos (solid curves in Fig. [5] differ from unity). 
This remains true, although the dependence on a is re- 
duced, if the models are normalized to have the same (a- 
dependent) linear theory rms fluctuations today (com- 
pare solid and dashed curves in Fig. [5]). If normalized 
to have the same (linear theory) rms fluctuations today, 
whatever the value of a, then the trends can be reversed 
(compare dotted curves with unity in Fig. [5]) . This last 
normalization convention is sometimes known as 'cluster- 
normalized': our work suggests that, in the context of 
modified gravity models, this jargon is misleading! 

We showed that our analysis captures the essence of 
the trends seen in the simulations (Figure [5]), so, in prin- 
ciple, the abundance of rich clusters should place inter- 
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esting constraints on modifications to the gravitational 
potential. In particular, the modification to cluster abun- 
dances cannot be reproduced by standard gravity with 
initial conditions modified to match the change in the lin- 
ear theory power spectrum; the differences in abundance 
can be larger than ten percent for sufficiently massive 
halos (Figure [5]) . However, to use cluster counts quanti- 
tatively in this way, our analysis should be extended to 
models in which objects form from an ellipsoidal collapse, 
as was necessary for standard gravity [H, [53] • 

Our analysis also helps to understand an interesting 
fact about the shape of the nonlinear power spectrum in 
modified gravity theories. Figure 7 in shows that 
whereas the large scale power spectrum in modified the- 
ories may be rather different than in standard gravity 
(because the linear growth factor is modified), the power 
on small scales (k > 1) is unchanged. Our analysis shows 
that, because small mass objects were never larger than 
the scale r 8 on which gravity was modified, they are not 
affected by the modification, so their abundance is the 
same as in the standard case. This is not affected by the 
initial conditions that we choose so that in both when a 
is non-zero, and when the power spectrum is changed so 
that we end up with a final power spectrum the same as 
in the case of modified gravity, the abundances of small 
halos is unaffected. In addition, their formation histo- 
ries will also be unchanged, so their internal structural 
parameters (shapes, density profiles) are also unchanged. 
In the halo model description of large scale structure [191 ] , 
the power at k > 1 is dominated by small mass halos. 
Since these are the ones for which gravity is essentially 
unmodified, the small-scale power spectrum is also un- 
changed. 

Figure [5] shows the result of a complementary study. 
In this case, we selected a massive halo from one of the 
simulations (say a — 0.5), and then looked at where its 
particles were in the other runs (with different a). The 
figure illustrates clearly that when a = 0.5, then the 
particle distribution is more compact. E.g., in the top 
panel, the large halo in the a = 0.5 run is broken up 
into three smaller haloes in the a = — 1 run. The bottom 
panel shows another effect: that the particles which made 



up a halo in the a = — 1 run are in a different location in 
the a = 0.5 run, suggesting that the peculiar velocities 
of the most massive halos are also sensitive to a. 

The sequence of contours associated with the different 
a runs look rather similar to the time evolution of an 
object in, say, an a = model. Thus, the most massive 
halos in models with a < may be like high-redshift ver- 
sions of the most massive halos in models with a > 0. 
Therefore, these figures suggest that the galaxy popula- 
tions in massive clusters may be rather different in mod- 
els with large a than when a is small. In particular, it 
is likely that the difference between the cluster and field 
galaxy populations increases as the strength of gravity 
on large scales increases. This is largely a consequence 
of the different values in these runs - so cluster M/L 
ratios, currently used to to constrain as [e.g. [5^, [3(| , may 
one day be used to constrain modifications to gravity. 

Furthermore, in standard gravity models, there is a 
strong correlation between evolution and environment 
[20j, |2l|, |28j . There are two reasons to suspect that this 
will be different if the gravitational potential is modi- 
fied. First, in the standard model, the correlation be- 
tween local environment and evolution is a consequence 
of Birkhoff 's Theorem. Birkhoff 's Theorem is lost when 
the force law is modified (it is this which modified our 
spherical evolution model from the standard case). And 
second, Figure [8] shows that the time scale for the as- 
sembly of objects is modified. The environmental de- 
pendence of galaxy properties is in rather good agree- 
ment with the standard model [3l|, HH , so it may be 
that this will one day provide interesting constraints on 
a. This is the subject of work in progress. In standard 
gravity, the formation and abundance of voids can be es- 
timated using similar methods to those used for clusters 
[34j - extending our analysis of the modified potential to 
voids is also work in progress. 
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APPENDIX A: FITTING FORMULAE FOR THE 
'MOVING' BARRIERS 

This Appendix provides fits to the barriers in Figure 
[31 First, let x — log M/(pr^), where M is measured in 
units of 10 10 /i _1 M & and r s is measured in /i _1 Mpc, then 



for positive a, 

55 ^ 5aJ -l-(^«U4)(Uf/* n a«'. (Al) 
whereas for negative a, 

5,5^-1 + |o^B|"(»ilf/* n fl-. (A2, 

These fits are accurate to a few percent in the range of 
(—7,3) for x, r s up to 20h~ 1 Mpc, and a between — 1 
and 1. For r s = 5/i _1 Mpc, the high end of the range 
of x corresponds to 10 16 /i _1 M Q . The dependence here 
makes sense, because we expect the mass scale of the 
modification to scale as r%. As for the a dependence, we 
can see from figure [3J that for negative a the deviation 
from standard gravity is stronger, hence the dependence 
on a for positive a, and j ck | 1 " 5 for negative a. 

In principle, one can use these expressions to gener- 
ate analytic approximations to the halo mass function 
as follows. For a given initial power spectrum, x can be 
written as a function of Si{M)\ this specifies the barrier 
shape in the units which are useful for the excursion set 
approach. Then insert this barrier shape into the expres- 
sions for the first crossing distribution given by (2l| . The 
ratio of this first crossing distribution to that associated 
with a barrier of height S c quantifies the change in halo 
abundances which is due to a and r s . Multiplying this ra- 
tio by the actual halo abundances in the standard model 
[2(| yields an analytic expression for the abundances in 
the modified model. 



